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Abstract 

Resistance to chemotherapies, particularly to anticancer treatments, is an increasing medical 
concern. Among the many mechanisms at work in cancers, one of the most important is the 
selection of ttimor cells expressing resistance genes or phenotypes. Motivated by the theory of 
mutation-selection in adaptive evolution, we propose a model based on a continuous variable that 
represents the expression level of a resistance gene (or genes, yielding a phenotype) influencing in 
healthy and tumor cells birth/death rates, effects of chemotherapies (both cytotoxic and cytostatic) 
and mutations. We extend previous work by demonstrating how qualitatively different actions of 
chemotherapeutic and cytostatic treatments may induce different levels of resistance. 

The mathematical interest of our study is in the formalism of constrained Hamilton-Jacobi 
equations in the framework of viscosity solutions. We derive the long-term temporal dynamics of 
the fittest traits in the regime of small mutations. In the context of adaptive cancer management, 
we also analyse whether an optimal drug level is better than the maximal tolerated dose. 

Key words: Mathematical oncology; Adaptive evolution; Hamilton-Jacobi equations; Integro-difFerential 
equations; Cancer; Drug resistance 
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1 Introduction 

Many pharmacotherapeutic processes, in anticancer, antiviral, antimalarial, or antibiotic therapies 
may fail to control proliferation, because the target (virus, cell, parasite) population becomes resistant 
to the drug(s). This may occur either because of evolution independent of drug(s), or because the 
drug(s) select for resistance. The occurrence of drug resistance is thus a major obstacle to therapy 

success. 
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Resistance is commonly seen in a diverse range of systems including agriculture pest resistance to 
toxic crops j3H I45j. bacterial resistance to antibiotics [1], mosquito resistance to antimalarial treat- 
ments [21 and references therein) and bacterial resistance to high copper concentrations [27j . 

There is overwhelming evidence that populations of cancer cells may harbour a diverse array of 
chemo- resistant subpopulations. 
Resistance mechanisms in cancer may be classified according to at least two major effects: 

• Cell adaptation. Cells may adapt metabolically in response to drugs. A typical example is the release 
of resistance proteins such as P-gp and other ABC transporters (ATP Binding Cassette), and studies 
on large human cohorts have shown their correlation with resistance in Acute Myeloid Leukemia [35] ) , 
since they promote the efflux of toxic molecules from leukaemic cells |TOl I24j . Induction of ABC trans- 
porter gene transcription by stressful conditions (including exposure to anticancer drugs) is a likely 
contributing mechanism, as reviewed in [IT]. This resistance mechanism can further be enhanced 
by protein exchange between cells [30], and its properties have been analysed using a mathematical 
model |3|ij. In as much as these mechanisms occur in individual cells and do not involve differential 
selection on variants resulting from mutations, we distinguish such cell adaptation mechanisms from 
other mechanisms stemming from mutations. 

• Mutations. It is generally accepted that tumor cells usually have higher division and higher mutation 
rates than normal cells; possible mechanisms for this have been investigated in bacterial populations 
in stressful conditions, inducing the formation of error-prone polymerase |30j . from which one may 
infer similar explanations for cancer cells. Moreover, when a drug targets a specific DNA site or a 
specific protein, a mutation that deletes this site or induces a change in protein conformation may 
confer resistance to mutant cells 1421. 



To elicit the relative contributions of these two mechanisms presents a considerable challenge. The 
answer may depend on many parameters such as the type of tumor and therapy. In the latter case, 
such mutations (possibly by a single amino acid of the target protein) may induce resistance to the 
drug at stake, whereas in the former case of multiply targeted drugs (or drugs with unidentied specie 
targets) a re- wiring of intracellular signaling networks is a plausible mechanism, as discussed in [33]. 
As a first approach to this complex question, we focus in this paper on continuous mechanisms that 
involve resistance induced by the drug environment on cell populations, and not on mechanisms that 
can be represented by a constant probability of mutation resulting in drug resistance, as in classical 
studies such as [22] or more recent work based on evolutionary models in a stochastic setting [29]. 

In this paper we are particularly interested in the latter mechanism, which is closely related to the 
process of Darwinian evolution, according to which certain rare mutants may have positive growth 
rates and be selected in environments that would otherwise result in ecological extinction. The term 
of 'evolutionary rescue' has been coined for situations where a population is saved from extinction 
by natural selection of genetic variants either existing in a population when environmental condi- 
tions deteriorate, or emerging during the otherwise inexorable decline to extinction. In this view, 
because some mutant cells have a fitness advantage, they will come to dominate dynamics and the 
population as a whole. We propose a mathematical model that takes into account birth, death and 
mutation rates of healthy and cancer cells depending on the level of resistant phenotype expression. 
Different expression levels would be refiected by the cell content of proteins (e.g., ABC transporters) 
responsible for expelling toxic molecules. We assume this expression level to be a continuous vari- 
able, and as such amenable to provide a basis for structured population models designed according 
to a selection/mutation type of formalism previously developed using Partial Difi^erential Equations 
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The phenomenon of drug resistance is usuaUy described mathematically by considering two distinct 
types of cells, sensitive and resistant. A recent probabilistic approach |18j incorporated dose modula- 
tion to study the evolution of resistance under various dosing strategies. Unlike probabilistic models, 
we do not consider that a cell or a cell population is necessarily either totally sensitive or totally 
resistant to a given drug; rather we introduce a continuous physiological variable describing resistance 
between (completely sensitive population) and +00 (completely resistant population). Such a vari- 
able, which may be genetic or epigenetic [l2], indeed provides the basis for a physiologically structured 
model of resistance. 

Mathematical models [21j incorporating phenotypic and microenvironmental factors suggest that 
adaptive treatment strategies can improve long-term survival compared to a maximum dose strategy. 
This may be achieved at the cost of maintaining the persistence of sensitive tumor cells, which are 
fitter than the resistant cells in low drug pressure conditions (i.e., not trying to eradicate them, an 
apparently paradoxical view ('adaptive therapy') [201 121], most likely applicable to slowly developing 
tumors only). Conversely, in cases of sudden rapid proliferation (e.g., acute leukaemias), only the 
employment of maximum tolerable doses is thought to have chances of controlling or eradicating the 
cancer cell population. 

Furthermore, depending on the nature of the tumor, different sorts of drugs can be used either 
separately or in combinations, in order to reduce the probability of resistance emergence and reduce 
side-effects on healthy cells |43| . This aspect is studied in [46^ W7\ using ordinary differential systems 
which give a macroscopic understanding on the use of multitherapies. In particular, besides cytotoxic 
drugs, here we focus on an additional class of therapeutic agents, so-called cytostatic drugs, which 
act by slowing down cancer cell proliferation and tumor growth. Cytostatic drugs have lower toxicity 
for healthy cells and reduce the emergence of resistance, which usually follows from treatments with 
cytotoxic drugs. In fact, they allow the survival of a small number of chemosensitive cells, which can 
reduce the growth of resistant clones through competition for space and resources. Moreover, whereas 
cytotoxic drugs remove sensitive cells from the population, they do not prevent the proliferation of 
chemo-resistant cells. There are two important effects of such therapies on cancer resurgence and 
its aggressiveness. Firstly, proliferating resistant cells may produce a series of increasingly aggressive 
subpopulations through natural selection p]. Secondly, toxic effects of the drug on chemoresistant 
cells, although not lethal, may induce genetic instabilities resulting in the generation of more aggres- 
sive mutant cancer cell lines [19]. These concerns within an expanding chemoresistant population 
are reduced under combined chemostatic and chemotoxic therapies, since toxicity of the latter is re- 
duced, and the chemostatic drug holds in check the multiplication of the chemoresistant subpopulation. 

The paper is organized as follows. We first present a model for healthy and tumor cells, structured 
by a variable representing the expression level of the resistance gene. We consider cases with and 
without chemotherapy. We also analyze the global behavior of solutions without mutations. Then, in 
Section [3j we turn to the analysis of situations with small mutations and employ the methodology of 
populational adaptive evolution to describe the dynamics of the fittest cells. Section |4] is devoted to 
the question of analyzing, for this simple model, whether the maximal tolerated dose is the optimal 
choice to impede the evolution of resistance. Section [5] presents several numerical simulations that 
illustrate our theory. Finally, we propose in Section |6] a model for the dynamics of healthy and cancer 
cells exposed to cytotoxic and cytostatic drugs. The model relies on the same structured population 
formalism presented at the beginning of the paper and it is analyzed by means of numerical simulations 
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looking for combined doses resulting in treatment optimization. 



2 A structured population model for healthy and tumor cells under 
the effects of cytotoxic drugs 

A major constraint in chemotherapy is to keep the toxic effect on healthy cells below a critical thresh- 
old. For that reason, it is preferable for a model to be able to represent both healthy and tumor 
cells. This is what we propose here in a simple equation that can take into account both types of 
cells, thus tackling the two major issues of medical treatments in general: adverse toxic effects in 
healthy cells and occurrence of drug resistance in diseased cells. At this stage we do not represent the 
spatial extension of tissues; this is an ingredient that can be included in a second stage. It is indeed 
more important from a therapeutic point of view to structure the two populations according to a 
physiological variable, e.g., a continuously evolving genetic trait or a cell content in specific proteins, 
that physiologically describes the evolution of the population as influenced by the treatment, than to 
include a space variable in the model. 



2.1 Selection/mutation model 

Let nH{x,t) I nc{x,t) denote the population density of healthy / cancer cells with gene resistance 
expression level x at time t. In the sequel we will use the term 'gene expression' meaning not only 
expression of one supposed resistance gene, but more generally of several genes yielding together a 
continuous drug resistance phenotype, more or less in the same way as others represent evolution to- 
wards malignancy in colorectal cancer [15], revisited by [H]. Note that in the model we develop below, 
we do not represent evolution towards malignancy, assuming the existence of an already constituted 
cancer cell population (with subscript C), opposed to a healthy cell population (with subscript H). 
The resistance level x to the drug can be measured either by the average molecular cell concentration 
or activity of ABC transporters that are known to be associated with resistance to the drug, or by 
the minimal drug concentration to kill the population (more precisely the LD50, LD90... of pharma- 
cologists, i.e., the minimal dose to kill a given percentage of the cell population). Different strains of 
the same initial cell lineage, selected by exposure to progressively increasing doses of the same drug, 
are indeed available in tumor banks, such as acute myeloblastic leukemia (AML) cells at the Tumor 
Bank of St. Antoine Hospital in Paris |48| . We a priori take < x < oo, even though the design of 
the model should lead to limitations on the possible values of x. 

The growth dynamics of healthy and tumor cells with a chemotherapy is given by the system 



growth with homeostasis 



natural apoptosis effect of drug 







(1) 



birth with mutation 
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dt 



nc{x, t) = (1 - 9c) r{x) - d{x) - c{t)nc{x) nc{x, t) 



(2) 



+0c J r{y)M„^{y,x)nc{y,t)dy, 
and the total population is defined as 

poo POO 

p{t) = pH{t) + pc{t), pH{t)= nH{x,t)dx, pc{t) = nc{x,t)dx. 

Jx=0 Jx=0 



(3) 



The following notations, interpretations and assumptions are used: 

• r{x) denotes the basic reproduction rate and d{x) > denotes the basic death rate, which depend 
on the gene expression level x. In order to incorporate a cost to produce the resistance gene, we 
assume r is decreasing, d is increasing 

r(0) > d{0) > 0, r'(-) < 0, r(+oo) = 0, d'{-) > 0. (4) 

• < Oh,c < 1 denotes the proportion of divisions with mutations and we can assume it to be higher 
for cancer cells. 

• /3 > is introduced, with the simplest possible form, to impose healthy tissue homeostasis (see 
below). For tumor cells, uncontrolled proliferation is obtained by birth terms that do not depend on 
the total population (i.e., we do not consider density-dependent inhibition in cancer cells). 

• c{t) denotes the dose of chemotherapy. Here we assume it has only an effect on increasing apoptosis. 

• tJ'H,c{x) represents the phenotypically dependent response to the drug; drugs are designed to target 
cancer cells more than healthy cells. The effect of drugs is here assumed to be summed up directly 
on mortality (i.e., in this simple setting, not involving the cell division cycle, we do not consider drug 
effects at cell cycle phase transitions in proliferating cell populations, see e.g., [28] for a discussion 
on this point) with a rate p that depends on the gene expression x (sensitivity to the drug) and we 
assume the properties 

IJ-H,c{-) > 0, '^^^''^ < 0, pH,c{oo) = 0, PH < PC- (5) 

This means that the drug always decreases the growth rate but its efficiency is lower on cells expressing 
a higher resistance gene expression. Also because the drug is efficient against cancer cells we assume 
that c is such that 

r(0) -d(0) -c^c(O) < 0. (6) 



• Mfjjj p (y, x) > denotes the probability that a mutation of a cell with gene expression y leads to 
a daughter cell with level x and (Th,c > measures the average size of these mutations. This means 
that 

fOO 

/ xMfj{y,x)dx = a. 

Jx=0 

We can assume that this size is larger in tumor than in healthy cells, in which control of the genome 
is more strictly ensured and mutations, if they occur, are of little consequence; hence ac > cth- To 
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be consistent with this interpretation, we impose (further conditions are imposed later on) 



poo 

/ M,{y,x)dx = l. (7) 

Jx=0 



It is useful to write the equations on nc and uh as 
d 



g^n{x,t) = n{x,t)R{x,t) +e J r{y)M„{y,x)[n{y,t) - n{x,t)]dy, 

so as to better see that the net growth rate (fitness) of the cells is given by 

RH{x,t) =- — r{x) - d{x) - cfiH{x), 

i^ + PiW (8) 
Rc{x,t) = r{x) — d{x) — cfic{x)- 

2.2 Healthy or tumor cells without mutations 

In order to motivate the introduction of these different terms, we will prove in the next sections that 
the model solutions exhibit behaviors that are in accordance with the biomedical interpretation. It 
is simpler to state results when there are no mutations, which we do now and postpone to the next 
sections the case with mutations. 

We introduce our assumptions and notations in each one of these three cases of interest indepen- 
dently: healthy cells only, cancer cells without therapy, resistance in cancer cells generated by therapy. 

Healthy tissue, no therapy, no tumor cells. In a healthy tissue we take nc = 0. Then, 
in the absence of therapy (i.e., c{t) = 0), the classical theory, that we recall below, shows that as a 
vanishes, the population model always selects the gene expression x = 0. In other words, x = is a 
globally attracting Evolutionary Stable Distribution: 

Lemma 2.1 (Healthy cells only, no mutations) We assume Q and that all data are Lipschitz 
continuous; we also assume 

eH = 0, n%>0, n%eL\ nc = 0. 

Then, we have 

nH{x,t) — y PH.ooKx), 

weakly in the sense of measures where the number density of cells ph,oo defined by zero growth rate 
(homeostasis) is computed as 

^ r(0) = d(0). (9) 



+ PH,oo) 



/9 



In words, with no mutations to create variability, the monoclonal population is globally attractive 
with the fittest trait x = 0. 

The proof of this result follows from a general analysis that can be found in \32\ 1371 EH] ■ It uses an 
estimate on the total variation of pnit) that holds in a very general framework that includes the case 
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at hand. 



Keeping mutations, a similar conclusion can be reached but in the framework of asymptotic analysis; 



to show this is a job that uses more elaborate ingredients and that is explained in Section 3.1 



Cancer cells, no therapy. For healthy cells, the total cell number p{t) induces a limitation 
on the multiplication of cells that represents homeostasis (i.e., the capacity of a tissue to stop growing, 
then differentiate and commit itself to fulfill a given physiological task). For cancer cells, we observe 
an uncontrolled growth in the absence of therapy whatever the gene expression level (at least close 
enough to x ~ 0, otherwise the cost to produce it can compensate the birth/death ratio). However, 
they remain with the lowest resistance gene because of the reproduction advantage when a; = 0. This 
can be derived from the model since we have the 

Lemma 2.2 (Cancer cells only, no mutations) When c{t) = and 6c vanishes, solutions to ^ 
with an initial data (x) > near x w satisfy, with an exponential rate 

Pc{t) — > oo, and weakly in the sense of measures '~^^,\^ — > 5{x). 

t—>-co Pc{t) *->oo 



Proof. From assumption Q, there is a value Xc > for which 

r(x) - d{x) > gm := ^[r(0) - d{0)] > 0, Vx G [0,Xe]. 
We write after integration in x 

nc{x,t)dx = / [r{x) — d{x)]nc{x,t)dx > gm / nc{x,t)dx. 
Jo Jo 

Therefore the first statement is proved since 

nc{x,t)dx > e^'"* / nc{x)dx, Xm = 9m- 



d 
~dt 



We also notice that the same argument, using a small interval closer to x ~ tells us that for all 
A < r(0) — d(0), there is a x\ such that 

nc{x,t)dx > e^^ / n^(x)(ix, 



and thus, for all A < r(0) — ti(0), there is a value p{\) such that 

Pc{t) > p{X)e^\ 

From this, we easily conclude the second statement which is a consequence of Laplace formula. Indeed, 
for a given y, we can choose from assumption (j4]), A > r(y) — d{y) and thus for x > y 

nc{x,t) , ei^(^)-^(^)^^ 
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In other words, the cell density is asymptotically small for all x > 0, and thus the solution concentrates 
at a; = 0. □ 
The case with mutations is treated in Section! 



Cancer cells with therapy This is the most interesting situation where resistance may follow from 
selection of cells with a high level of gene expression. This can arise because the drug concentration c 
is limited in order to keep healthy cells below a given toxicity threshold. This situation occurs if the 
maximum net growth rate of tumor cells, as defined through ([s]), is positive for some x > 0. Here we 
concentrate on the simplest case of a constant therapy and we make the additional assumption that, 
for all level of c, there is a maximal resistance gene expression Xc > which maximizes the fitness 

for X 7^ Xc, r{x) — d{x) — cfic{x) < Rc '■= r{xc) — d{xc) — cfic{xc) > 0. (10) 

Notice however that 

R~c < r{xc) - d{xc) < r(0) - d{0) = Rq. 
With this assumption, we can establish the 



Lemma 2.3 (Cancer cells with therapy, no mutations) Assume c = I, 6c = and (10). Then, 
the solutions to ^ with an initial data n^{x) > satisfy, with an exponential rate 

f+\ . nc{x,t) 
Pc{t) ^ — ^ oo, -— — y6{x-Xc). 



We do not repeat the proof of this result which is the same as for Lemma 2.2 and gives that 
pc{'t)e^^ — > oo for all A < Rc- 



t—^ca 



3 Asymptotic analysis with small mutations 

A clear mathematical way to precisely express the results mentioned above, i.e., that a specific level 
of gene expression is selected, is through asymptotic analysis for small mutations and observing the 
dynamics in the long run. In ecology, these models generate populations that are highly concentrated 
on some well separated traits with possible branching and polymorphism. For the case at hand, the 
mutation rates might be too high to reduce the analysis to this limit, but it is still a way to express 
mathematically what happens for small, but non zero, frequencies of mutation. 
In this section we state the behaviors that can be derived from analysis in the different cases already 



mentioned without mutations in Section 2.2, They give details on the evolution of the population 
towards the ESD (Evolutionary Stable Distribution, see |26j). which is that at each time a distribution 
close to a Dirac mass is obtained. 

Even though these are not the most interesting cases, the results listed below are also valid when 
there are no mutations (Oh = or Oq = 0) and this is the case used later in the numerical simulations 
of the solutions. 

3.1 Asymptotic analysis: case of healthy cells 

With the aim of studying the dynamics of the model in the aforementioned limit of small mutations 
and many generations, we introduce a small parameter e, which is used to rescale time and define a 
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more specific choice of the mutation kernel, so that the population dynamics ([T]) for healthy cells only 
and without treatment can be rewritten in the form 



d 

e—nH{x,t) 



I -Oh 



[l + PHit])' 



r{x) — d{x) nnix, t) 



+ - 



roo 

PH{t)= / nH{x,t)dx. 
Jx=0 

We have written the mutation kernel in such a way that 



6h /" / n 1 T7/ y ~ X-. , ■. , 

p r{y)-M{y, )nH{y,t)dy, (n) 



Ma„{y,x) = lM{y,^) with sup / e^^^"" M{x, z)dz < oo, 

x>0 . 



(12) 



M{x, z)dz = 1, iZ'^ zM{x, z)dz = 0, 



so that it appears directly that x — y = 0(e), in other words, mutations have a small effect in terms of 
the change in the x, y variables. Notice that, in these expressions, we always have ez < x. Therefore, 
the last assumption implies that cells in the state y = do not undergo mutations (they could only 
be to z < thus contradicting the zero integral condition). This is useful to state clear and simple 
results, otherwise the evolutionary drift creates an ESS (Evolutionary Stable Strategy) (see [13] for 
an introduction to the dynamical system theory of adaptive dynamics) for a slightly positive x and 
this can be analyzed mathematically with our method but we will not do it here for simplicity of the 
presentation. 

The analysis in [UISIEH] can be used and proves that, thanks to the assumptions Q, the fittest pop- 
ulation is obtained when the gene is not expressed at all, in mathematical words, nnix, t) = ph,oo^{x) 
is attractive for this system. 

Before stating a precise result, we give assumptions. Firstly, we assume that the initial population is 
concentrated as a 'sharp Gaussian' around a state corresponding to a non-vanishing gene expression xf' 



nH,e{x,0) = e^H.e'y^)!'^ ^ m^^(x) — > u^ni^) locally uniformly 



£->0 



max'u^(x) = = u^jj{x^), nH,e{x, 0) — > P%^{x — x^). 



(13) 



e->0 



Notice that, as e vanishes, only mutations can deviate such a population determined by p^6{x — x^) 
from this initial state which is asymptotically a steady state for 9h = as soon as the total population 
satisfies 

^ r(x°) = (i(x°). (14) 

(1 + pO^)^ 

More precisely, with these new notations we may write 



Theorem 3.1 (Healthy cells only) We assume Q, (12)-(14) and that all data are Lipschitz con- 
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tinuous. Then, the population is monoclonal in the limit of small mutations: 

nH,eix, t) — > pH{t)5{x - x{t)), 

e-s>0 

and it returns to the fittest trait x = 0: using the notation ([9]), we have 

x{t) — > 0, pnit) ^ — > Ph,oo- 



Proof. The limit e — t- 0. The proof is based on the use of the Hopf-Cole transform 

duf, e duH 



Ue{x,t) = eln{nH{x,t)), 



dt uh dt 



Insertmg it in (11), we find the equation 



= / rix) - d{x) + ^ . / r{y)M{y, ^)e("(^.*)-(-.*))/^dy 

and changing the variable y ^ z with y = x + ez, we find 

= „ r{x) - d{x) + — ^ / r{x + ez)M{x + ez, ^)e(«-(^+=^.*)-«-(^'*))Adz. 

At this stage, it is useful to introduce the limit in the integral term which is given by the Hamiltonian 

H{x,p) = r{x) [ M{x,z)eP-'dz. (15) 



Notice that from (12), we derive the properties 

H{x,0) = r{x), Hp{x,0) = 0, Hpp{x,p)>0. 
By convexity, these imply in particular 

H{x,p) > r(x). (16) 

One can establish two types of uniform estimates: (i) the total number of cells pH,e{i) is uniformly 
bounded in BV{Q^T) for all T > and {■^PH,e{t)) _ is bounded, therefore 

PHAt) PHit), j^PHit) > -c. (17) 

(ii) Uniform Lipschitz estimates can be established for this kind of models (see [Ij for the specific case 
of integral mutations). 

Therefore, we may now pass to the limit e — )• in the equation on Us in the viscosity sense (see 
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references below), and we find 



d 



I -Oh 



dt {l + pH{t)y 
max u{x, t) = 

0<x<oo 



r{x) — d{x) + 



Vt > 0. 



Oh 



H{x,\7u{x,t)), 



(18) 



This equation is the fundamental new object that describes the dynamics for e — )• 0. It is a constrained 
Hamilton-Jacobi equation, which is rather different from the usual Hamilton-Jacobi equation because 
the L°° contraction property is lost for instance. Nevertheless, the correct notion of solutions are still 
the viscosity solutions of Crandall-Lions (see [3l IIH [T6]). 



The solution to equation (18) is a pair {u,p) where p{t) is the Lagrange multiplier associated with 
the constraints maxx u{x,t) = 0. There are two possible approaches to rigorously prove this limit. In 
the approach developed in [32], we consider strong assumptions which allow for uniform smoothness 
of Us, and then the theory can be carried much further including the long time dynamics of x{t). In 
this approach concavity assumptions are needed which are difficult to meet in the case at hand (it is 
an interesting problem to develop the theory). 

The other approach, that we follow here, has been developed in [U [5l [HI [38] . We consider a weak 



theory with our general assumptions and viscosity solutions to (18). The function u is merely Lipschitz 
continuous and semi-concave [38], the best possible regularity is p{t) £ BV{U.'^) because jumps on 
p{t) and x{t) are possible. This analysis for weak solutions (only possible in one dimension) gives also 
that, along with the dynamics, the fittest trait is characterized by 



= u{x{t),t), 







1 



r{x{t)) - d{x{t)). 



(19) 



This second equality is just a consequence of ll = Vu = at a maximum point (viscosity solutions 
are not enough and semi-concavity is necessary here). 



The limit t — )• oo. Additionally, we also know from (17) that pnit) is non-decreasing, therefore taking 
into account (19) and the definition ([9]), we find 



PH{t) pH < PH,i 

t—>oo 



X{t) \ 
t— >oo 



It remains to identify these limits. 
By contradiction, assume that Ph < Ph oo, then we use that 



Ft 



u{x, t) 



> 



1-Jh 

OH 
1 



-^r{x)-d{x) + 
^ + PHit)r 

r(x) — d{x) 



Oh 

[^ + PHit)) 



^ H{x,Vu{x,t)) 



'^+PHit)Y 



using (16) for the inequality. With our assumption and the definition of pH,oa, there are intervals 



[xi,a;2] where the right hand side is positive for all t and thus ^u{x,t) > (uniformly) on this 
interval which contradicts the condition max^ u{x, t) = 0. This proves that pn = Ph,oo and concludes 
the proof of the theorem. n 
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3.2 Asymptotic analysis: cancer cells, no therapy 



Following the lines drawn in Section 3.1 , we can arrive at the same conclusion that, considering cancer 
cells only, the gene expression x = is selected even when mutations are present but small. Again a 
clean mathematical path towards this direction is to rescale the model for cancer cells and, following 



the lines of Section 3.1, to rewrite it as 



(1 - Be) r{x) - d{x) ncix, t) + % 



)nc{y,t)dy. 



(20) 



Because it is linear and expresses unlimited growth, this model is very different from those in Section 



3.1 for healthy cells (or in ecology). In these cases there is a limitation by nutrients or space which 
controls nonlinearly the growth and the dynamics through an integral variable p{t). Nevertheless 
we can use the same method to analyze again the main effect, which is that x = is a globally 
attractive ESD; the proof highlights how a new control unknown shows up naturally and gives the 



same constraint and Lagrange multiplier in the limit as in (18). To state the result we need again to 



make more precise that the initial data is supposed monoclonal at a value xq > 



maxM^(x) = 



(x) 



[x) locally uniformly, 



n^(xO), 



nc eix, 0) — p^5(x — x^). 

' e-s>0 



(21) 



Theorem 3.2 (Cancer cells, no therapy) We assume Q, (12), (21) and that all data are Lips- 
chitz continuous. Then, the population is monoclonal in the limit of small mutations: 



ncix,t) 



/o°° nc{x, t)dx £->-o 



5{x — x{t)), 



and it returns to the fittest trait x = 



x{t) — y 0. 

t— >oo 



The model again is not very relevant for direct interpretation because, as we see it below, the cell 
number density can decrease exponentially and then it explodes exponentially. A behavior reflecting 
only the fact that cells cannot be naturally endowed with a gene expression x ^ 0; this may occur if 
a treatment is stopped after resistance has been triggered (see next section). 

Proof. We define the probability density 

, , nc(x,t) . . f°° . , , 

Pe[x,t) := 7-^, Pc{t) = nc{x,t)dx. 

Pc{t) 







Because pc{t) satisfies the equation 

^ dpc{t) 
^ dt 



oo 



[r(x) — (i(x)]nc(x, t)dx 







12 



we also find a closed equation for p^, namely 



d 



1 



e^p,(x, t) = peix, t)[{l - ec)r{x) - dix)] - pe{x, t)Ie{t) + Oc / r {y)p ,{y) -^M {y , ^-^)dy 



£ 



with 



/•oo 

,{t) := / peiy,t){riy) - d{y))dy. 
Jo 



(22) 



As mentioned earlier, this integral comes naturally to play the role of the population density p{t) in 
the case of healthy cells as a Lagrange multiplier and we are back to the general setting in [H [5l [32l [38] . 
We can perform the Hopf-Cole transform again Ue = elnp^ and write 

^n, = (1 - 0H)r{x) - d{x) - Is{t) + eH I r{x + £z)M{x + ez, ^)e(«-(^+=-.*)-«-(^''*))/=dz. 



We can finally pass to the limit and find, still with the Hamiltonian defined in (15), the constrained 
Hamilton-Jacobi equation 



r d_ 

dt 



u = {l- eH)r{x) - d{x) - I{t) + Oh H{x, Vu{x, t)), 



max u{x, t) = Vt > 0. 

0<x<oo 



(23) 



As mentioned in Section [3. 1[ we find again that the fittest trait x{t) is characterized by 

= u{x{t),t), = r{x{t)) - d{x{t)) - I{t). 

From the definition of leit), this last equality is of course equivalent to the limit Pe — )• 6{x — x{t)). 

Again, the general theory applies here and we may conclude the proof as in Section 3.1 I{t) is 
non-decreasing and thus x{t) is non-increasing; this proves that these two quantities have limits for 
t — >• oo, which we easily identify from the Hamilton-Jacobi equation. 

With this information, notice that the dynamics of pcit) can be approximated by 

dpcit) 



dt 



pcitMxit)) - dixit))], 



and it may decrease exponentially if rix^) — dix^) > 0, but eventually it will increase because we know 
that r(0) - d(0) < 0. □ 



3.3 Asymptotic analysis: cancer cells with therapy 

In the case when treatment is included, the dynamics of cancer cells given by equation ([2])-(|3]) is 
independent from healthy cells and we concentrate on them. We use the same notations as before 
to express the small mutations regime where we can state a clear mathematical result and write the 
dynamics as 



£i^ncix,t) 



(1 - 9c) r(x) - dix) - pcix) ncix, t) + 9c 



r(y)^M(y,^ 



)nciy,t)dy. (24) 
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Theorem 3.3 (Cancer cells with therapy) We assume (10), (12), (21) forsomeO < xq < Xc 
and that all data are Lipschitz continuous. Then, the population is monoclonal in the limit of small 
mutations: 

nc{x,t) ^ -^^^^^ 

Jq nc[x, t)ax e-s>o 

and it selects asymptotically the fittest trait x = Xc 



x{t) — ^ Xc 
t^oa 



The structure of this problem is Unear as in Section 3.2 Therefore we have to adapt the general 
theory of concentration as we did it for cancer cells only. There is no special difficulty compared to 
the proof of Theorem |3 . 2| and we skip it. We just mention that, again, we cannot analyze directly the 
cell population density Uc, because it blows- up exponentially, and the natural quantity on which we 
can develop our method is the probability density 



oo 



, , nc(x, t) / \ I / \ , 

Pe{x,t) := 7—, Pc{t) = nc{x,t)dx. 

Pc{t) 







4 A consequence for optimal therapy 



As we have seen in Theorem 3.3, the model predicts selection of a resistance gene whose level of 
expression Xc depends on the drug dose. This leads to the question to know if cancer growth can be 
minimized by using an optimal dose c. At this stage we only consider the case of a constant dose c. 
We analyze the particular case of coefficients used for our numerical studies in Section [5] 

To state the problem, we recall that the resistant population will result in net growth (fitness) Rc 



given by formula (10) that is 



Rc := r{xc) - d{xc) - cudxc) = max [r{x) - d{x) - c^c(a^)] • (25) 

x>0 

Therapy optimization means finding c* that minimizes Rc over c > 0. The most usual protocol is to 
maximize the dose c with the only limitation of toxicity on healthy cells; does this minimize i?c? 
We analyze this question for the choice 

2 2 

Rdx):=-^-d-^-^; (26) 

here we have set := c to simplify the following calculations. Natural assumptions are 

0<a2(rg-d) <a^ (27) 
a < 1. (28) 

The first one is equivalent to Rc{x = 0) < 0, i.e., the therapy is high enough to kill the cancer cells 
in their normal state x = while cancer cells by themselves have a positive growth. The second one 
says that the therapy acts faster on the gene expression than the natural birth process. 
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In order to find the maximum of Rc{x), we use the variable y = x^, take the derivative and set it 
equal to zero 

K{y) = - , "\, + , f ., =^^ — = ^- 

^^^^ {l+yY (a2+y)2 l + y a?+y 

From this we obtain the condition for existence of a maximum of Rc{y) 

a^ro < a <^ R'^{0) > 0, (29) 

and the following expression for the maximum point 



Uc = = s 

^ otherwise. 

When yc > 0, we compute the maximum value of Rc 



(30) 



i^=^^P#-rf. (31) 
1 — 



We can now conclude that: 



• when a > ro, then Rc{y) increases from i?c(0) < —d to Rc = —d; this level of therapy is strong 
enough and resistance cannot occur. 

• When < a'^{rQ — d) < a < aPrQ (this interval of a is not empty for ro < 1 for instance), then Rc{y) 
decreases from Rc = Rc{0) > —d to —d. Although weaker, this level of therapy avoids resistance to 
occur. 

• When a^ro < a < ro, then Rc{y) increases from i?c(0) < to Rc and then decreasing; Rc can be 
positive iff (1 — a^)rQ > d and this can happen with Rc{0) < only when a^ro < ay^rg — d. Resistance 
occurs therefore if for a = a\J — d, then Rc > 0. An example is explicitly computed in Section [sj 



However, since Rc in (31) is decreasing in a in the range under consideration here, the only way to 



avoid resistance is to increase the dose. 



5 Numerics 



In this section we illustrate the analytic results obtained before by numerical simulations performed 
in Matlab. For all these computations we use 4000 points on the interval [0, 1]. Moreover, we treat 
the case without mutations, i.e.. Oh = Oc = 0. The parameter e = 10"^ is used for the initial data 



which we take under the 'highly concentrated form' ( 13 ), (21 ) so as to obtain very localized solutions 
ncH according to the theory in |32j . 



We first illustrate the case of healthy cells. Figure 5.1 shows the time dynamics of the concentration 



point in (11) with 

r(x):=^, d:=0.4 

and the initial data 

nQ{x) := Co exp(-(x - O.lf/e), 
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where the constant Co is adjusted to enforce ||'i-o||l1([o,i]) = 1- The computations were done using 
a imphcit-exphcit finite difference scheme. It can be seen that it takes very short time to adapt to 
R{x(t), p{t)) = 0, that n{t,x) concentrates and its maximum moves towards 0. 



0.9 - 




Figure 5.1: (Healthy cells) The plot to the left shows the level-sets of n{t,x). The dotted line shows 
the unique solution x to R{x,p{t)) = as a function of t. On the right, we display the values of 
nH{t,x) at the end of the computation. Notice that nnit^x) concentrates and its maximum moves 
closer towards for larger time. 



Next we illustrate the case when the therapy induces resistance. Figure 5.2 shows the time dynamics 



of the concentration point in (24) with 

r{x) :- 



d:= 0.245, fi{x) :-- 



0.55^ 



0.52 + a;2 



1+X2' 

and initial data as before but centered at = 0.5 

710(2;) := Co exp(— (x — 0.5)^/e). 
As shown in Section HI the concentration point for large time can be computed explicitly: 



xc 



a — 



0.8165. 



1 — a 

Since we have exponential growth of the maximum of ric, after a certain time the computations are 

breaking down. In or( 
,. n, 
equation on 



and Figure 
defined by (l22])). 



/q°° He dx 



are shown in Fig 


5ure 


5.3 


Figure 


5.2 


- d{x) - cucit) - 


m 




(with 1 
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Figure 5.2: (Resistance) The plot to the left shows the level-sets of nc(t, x)/ j nc{t, x) dx. The smaller 
of the two positive solutions of r{x) — d{x) — cucit) — I{t) = is shown as dotted line. On the right, 
we plot the values of ncit,x) at the end of the computation. Calculations are done for the equation 



(24) on nc- 




0.2 0.3 



0.8 0.9 



Figure 5.3: (Resistance ct'd) The plot to the left shows the level-sets of p. The smaller of the two 
positive solutions of r{x) — d{x) — cficit) — I{t) = is shown as dotted line. On the right, we plot the 
values of p{t, x) at the end of the computation. Calculations are done for the equation on p^. 



6 A structured population model for healthy and tumor cells under 
the effects of cytotoxic and cytostatic drugs 

We turn to the study of the coupled dynamics of healthy and tumor cells exposed to both cytotoxic and 
cytostatic drugs. Our model relies on the same structured population formalism presented in Section 



2.1, However, to allow a better qualitative agreement with biological reality, some modifications are 
introduced, which make the model more difficult to handle analytically. Therefore, the behavior of its 
solutions is studied only through numerical simulations, whose results are discussed at the end of this 
section. 



As in Section 2.1, the number densities of healthy and cancer cells with resistant gene expression 
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level < X < oo at time t > are described by functions nH{t,x) > and nc{t,x) > 0. The 
selection/mutation dynamics of cells under the effects of cytotoxic and cytostatic drugs are given by 
the system below: 

mutations and renewal 



d 

nH{x,t) = — ( I rH{y)M{y,x)nH{t,y)dy - rH{x)nH{x,t) 



dt l + aHC2{t) 

1 + aHC2{t) 



+ [t^^^^^^^^ -dH{x)lH{t)]nHix,t) - ci{t)fiH{x)nH{x,t\ ,(32) 

effect of cytotoxic therapies 



growth with cytostatic therapies and death 



— nc(x,t) = — — ^ (I rciy)M{y,x)nc{t,y)dy - rc{x)ncix,t) 
at 1 + acC2{t) 

rcjx) 
1 + acC2{t) 



+ ( 1 ."^^^^^.N - dc{x)Ic{t) ] nc{x, t) - ci{t)ficix)nc{x, t), (33) 



where 



Init) ■■= aHHPH{t) + ancPcit), Ic{t) ■= acHPnit) + accPc{t). (34) 
Apart from integrals Ih and Ic, which have a different meaning with respect to the integral / consid- 



ered within the previous sections, the same notations of Section 2.1 are used in the above equations. 
In addition: 

• Kernel M{y, x) denotes the probability that a mutation of a healthy/cancer cell with gene expression 
y leads to a daughter cell with level x. 

• Functions rnix) and rc{x) stand, respectively, for the proliferation rate of healthy and cancer cells 
with gene expression x. Factors 

and 



1 + a_ffC2(t) l + acC2{t) 

mimic the effects of cytostatic drugs, whose concentration at time t is described by function C2{t). In 
fact, as previously noted, such therapies act by slowing down cell proliferation, rather than by killing 
cells. The average sensitivities of healthy and cancer cells to these drugs are modeled by parameters 
Oh, ac G K^- 

• Function dn^x) represents the death rate of healthy cells due both to apoptosis and deprivation of 
resources (e.g., oxygen and glucose) by cancer cells. Analogously, function dc{x) models the death 
rate of cancer cells because of the competition for space and resources with the other cells. Two main 



differences are here introduced with respect to the model presented in Section 2.1 , On the one side, we 
assume the growth of cancer cells to be hampered by the competition for space and resources among 
themselves and healthy cells as well. As a result, we multiply dc{x) by Ic{t)- On the other side, 
we consider apoptosis and competition phenomena involving normal cells as mediated by interactions 
with the surrounding cells. Therefore, instead of multiplying function rnix) by as in ([T]), we 

multiply dnix) by I nit)- It is worth noticing that both modeling strategies mimic the same effects of 
net proliferation and death. 

• Function ci(t) denotes the concentration, at time t, of cytotoxic agents, which are assumed to kill 
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both healthy and cancer cells with gene expression x at rates fJ-nix) and ncix), respectively. 

• Coefficients anc ^'iid uhh stand, respectively, for the average interaction rates between healthy and 

cancer cells or among healthy cells themselves. Analogous considerations hold for acH and acc- 

As in Section 2.1, model (32), (33) can be recast in the equivalent form given hereafter, in order to 
highlight the role played by the net growth rates of healthy and cancer cells, which are described by 
functional Rh{Ih,ci-,C2-,x) and Rc{Ic,ci,C2,x): 



d_ 
di 

d_ 
dt 



nH{x,t) = RHilHit),Ci{t),C2it),x)nH{x,t) + 

nc{x,t) = Rc{Ic{i),ci{t),C2{t),x)nc{x,t) + 



Oh 



1 + aHC2{t) 



with 



RH{lH{t),Ci{t),C2{t),x) 

Rc{Icit),ciit),C2it),x) 



1 + acC2{t) 
rH{x)il-9M) 



rH{y)M{y,x)nH{t,y)dy 
rciy)M{y,x)nc{t,y)dy, 



1 + aHC2{t) 

rc{x){l-9M) 
1 + acC2{t) 



dH{x)lH{t) - flH{x)ci{t), 

dc{x)Ic{t) - ixc{x)ci{t). 



The following considerations and hypothesis are assumed to hold: 
Parameters 9h,c and functions ^h,c satisfy the assumptions of Section 2.1 



• We consider intra-population interactions as occurring at a higher rate than inter-population ones. 
As a consequence, the interaction rates are assumed to be real numbers such that 



aHH, acc > 0, auc, acH > 0, 



o-HH > anc, acc > acH- 



(35) 



• In analogy with Section 2.1, with the aim of translating in mathematical terms the idea that 
producing resistance genes implies resource allocation both for healthy and cancer cells (eventhough 
this is debated), we assume functions rn and rc to be decreasing 

rH,c{-) > 0, r'{-)H,c < 0. (36) 

In order to mimic the fact that mutations conferring resistance to therapies may also provide cells 
with stronger competitive abilities, functions dn and dc are assumed to be decreasing 

dH,ci-) > 0, d'{-)H,c < 0. (37) 



• Since we assume that cytostatic agents are designed to be more effective against cancer cells rather 
against healthy cells, we make the additional assumption 

an < ac- (38) 

We analyse the model with these new ingredients through numerical simulations illustrating how the 
outputs can be influenced by different concentrations of cytotoxic and cytostatic agents. From a 
biological perspective, this means to use the present model as an in silico laboratory to highlight some 
mechanisms that may play a key role in the development of cancer resistance to therapies, with the 
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aim of providing support to the design of optimal therapeutic strategies. The following contents are 
organized into four subsections. Subsection |6.1| presents those assumptions that define the general 



setup for numerical simulations. Subsection 6.2 and subsection 6.3 are devoted to study the separate 



effects of cytotoxic and cytostatic drugs on the dynamics of healthy and cancer cells. Finally, subsection 



6.4 focuses on the combined action of these two classes of anti-cancer agents, looking for the existence 



of suitable doses allowing the design of optimal therapeutic strategies. 

6.1 Setup for numerical simulations 

Numerical simulations are performed as in Section [5] with 2000 points on the interval [0, 1]. Interval 
[0, T] with T = 2000dt is selected as time domain, where the unit time dt is chosen equal to 0.1. 
We choose the initial conditions 

unit = 0,x)= nc{t = 0, x) = n°(x) := C° exp(-(x - 0.5) Ve), (39) 

where is a positive real constant such that 

QH{t = 0) + Qc{t = 0) « 1. 

Parameter e is set equal to 0.01 to mimic a biological scenario where most of the cells are characterized 
by the resistant gene expression level corresponding to x = 0.5 at the beginning of observations. 
Assumptions and definitions given hereafter are used along all simulations: 

M(y,x) :=C7Mexp(-(y-x)V(0.01)2), Cm I exp(-(y - x) V(0.01)2)dx = 1, VyG[0,l], 

1.5 3 

9h = Oc = := 0.1, ruix) := ' „ , rc(x) := — — ^, a// := 0.01, ac-=l, 

1 + x^ 1 + x^ 

auH = acc = a := 1, anc ■= 0.07 acH ■= 0.01, 
dnix) := 0.5(1 - O.lx), dcix) := 0.5(1 - 0.3x), 

0.2 0.4 

^= (0.7)2 + x2' ^c{x) := (0.7)2 + ^2- 

Functions ci(t) and C2(i) are assumed to be constant, i.e. 

ci(t) := ci G M+, C2(t) := C2 G M+. 
The values of parameters ci and C2 are chosen case by case according to our aim in each subsection. 

6.2 Effects of cytotoxic drugs 

At first, we study the effects of cytotoxic agents only on the dynamics of healthy and cancer cells. 
With this aim, we run simulations for different values of parameter ci with C2 = 0. The obtained 



results are summarized by Figure 6.4 which depicts a scenario that is in good agreement with clinical 
observations. 

In absence of therapeutic agents (i.e., when ci = C2 = 0), we observe, at the end of the simulations, 
the selection for those cells that are characterized by a low expression level of the resistant genes 



and thus, due to assumption (36), by a strong proliferative potential. On the other hand, as long 
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as ci increases, the number of cancer cells inside the sample at the end of the observation window 
becomes smaller but, due to the typical side-effects of cytotoxic agents, even the number of healthy 



cells decreases. Furthermore, Figure 6.4 highlights how cytotoxic agents alone favor the selection of 
resistant clones. In fact, the average resistant gene expression level of cells tends to increase with the 
concentrations of therapeutic agents. 



n (t=2000,x) 



0^=0,0^=0 
_c^=1.75, 0^=0 
c^=3.5, 0^=0 




n (t=2000,x) 



■ / 

y 









=1.75,0^=0 




=3.5, 0^=0 




Figure 6.4: (Cytotoxic drugs only) Trends of nc{x, t) (left) and nn^x, t) (right) at t = 2000 for ci = 
(dashed-dotted lines), ci = 1.75 (dashed lines) and ci = 3.5 (solid lines). In all cases, parameter C2 is 
set equal to zero. As parameter ci increases, functions nc{t = 2000, x) and n//(t = 2000,3;) tend to 
be highly concentrated around some increasing values of x and their maximum values become smaller; 
this indicates higher resistance with higher doses of drug. 



6.3 Effects of cytostatic drugs 

Focusing on the action of cytostatic agents alone, we run simulations increasing parameter C2 and 



keeping ci = 0. The results presented in Figure 6.5 highlight the capability of the present model 



to mimic the effects typically induced by cytostatic agents on cells dynamics. In fact, under the 
considered values of C2, the dynamics of healthy cells is kept unaltered while the proliferation of 
cancer cells is reduced as long as the drug concentration increases. Furthermore, it should be noted 
that increasing values of C2 tend to slow down the dynamics of cancer cells, i.e., for larger values of 
C2, function nc{x,t) at the end of simulations stays closer to the initial data nc{t = 0,x). 

6.4 Combined effects of cytotoxic and cytostatic drugs 

Finally, we are interested in analyzing the effects on cell dynamics caused by the simultaneous action 
of cytotoxic and cytostatic drugs. As a result, we perform simulations with various values of both ci 
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n (t=2000,x) 



n (t=2000,x) 




c 


=0,c^=1 




=0,c^=3 




=0,c^=7 




=0, 0^=1 
=0, 0^=3 
=0, c„=7 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 6.5: (Cytostatic drugs only) Trends of nc{x,t) (left) and nH{x,t) (right) aX t = 2000 for 
C2 = 1 (dashed-dotted lines), C2 = 3 (dashed lines) and C2 = 7 (solid lines). In all cases, parameter 
ci is set equal to zero. Increasing values of parameter C2 lead the qualitative behavior of function 
nc{t = 2000, x) to become closer to the one of nc{t = 0, x). On the other hand, the trend of function 
nH{t = 2000, x) remains basically unaltered. 



and C2- The obtained results are summarized by Figure 6.6 and Figure 6.7 They show that there are 
values of ci and C2 leading to extinction of the cancer cells while keeping alive about one half of the 
healthy cells at the end of computations. This is consistent with experimental observations suggesting 
that different therapies in combination can avoid the emergence of resistance and minimize side-effects 
on healthy cells. Furthermore, it should be noted that such values of ci and C2 are smaller than the 
highest ones considered in previous simulations, which do not even allow a complete eradication of 
cancer cells (see Figure 6.4 and Figure 6.5). Therefore, these results also suggest that optimized 
anti-cancer treatments can be designed making use of proper combinations between cytotoxic and 
cytostatic agents, thus supporting the idea that looking for protocols based on treatment optimization 
can be a more effective strategy for fighting cancer rather than only using high drug doses. 



7 Conclusion and perspectives 

Based on a model adapted from ecology, we have presented a mathematical analysis of resistance 
mechanisms to drugs under the assumption that resistance is induced by adaptation to drug environ- 
mental pressures. We have used physiologically structured equations where the variable represents 
the gene expression levels. Theories developed in other contexts of population biology and Darwinian 
evolution underlay our analysis and finding that therapy can trigger emergence of a resistance gene 
and result in disease relapse. 
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nc(t,x) n^(t,x) 




Figure 6.6: (Cytotoxic and cytostatic drugs) Dynamics of nc{x,t) for ci = C2 = (top-left), ci = 
C2 = 1 (top-right), ci = C2 = 1.5 (bottom-left) and ci = C2 = 2 (bottom-right). As long as parameters 
ci and C2 increase, the maximum value of nc{t = 2000,2;) becomes smaller so that, under the choice 
ci = C2 = 2, function nc{x,t) tends to zero across time. 



Importantly, results suggest the feasibility of optimized anti-cancer treatments, based on specific 
combinations of cytotoxic and cytostatic drugs (see Section [6]). On the other hand, if only cytotoxic 
drugs are available, then (not surprisingly) the optimal therapy is always to administer the maximal 
tolerated dose (see Section [4]) . This scenario is somewhat different from that studied in |2T| , where 
the resistant and susceptible populations are in competition, a polymorphism that cannot occur in 
our model because it contains a single environmental variable (Cause competitive exclusion principle). 
It is worth noting that we have analyzed the possibility for optimization only in the restrictive case 
of a constant level of therapy and a possible future extension would be to modulate dose with time [T7] . 

The model is very simple at this stage and several directions for improvements are possible. As 
mentioned above, because there is only a single environmental unknown (the drug pressure and only 
one trait depending on it), polymorphism is not possible and thus a more realistic description of the 
micro-environment of the cell populations is needed (see [T0| and the references therein). Another 
research direction would be to include several qualitatively different drug resistance mechanisms, to- 
gether with realistic, targeted, and clinically acceptable multi-drug therapies. 
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Figure 6.7: (Cytotoxic and cytostatic drugs) Dynamics of nH{x,t) for ci = C2 = (top-left), ci = 
C2 = 1 (top-right), ci = C2 = 1.5 (bottom-left) and ci = C2 = 2 (bottom-right). As long as parameters 
ci and C2 increase, the maximum value of nnit = 2000, x) becomes smaller but about one half of the 
healthy cells is still alive at the end of computations. 
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A Matlab code for Figure 5.1 



°/o solves equation 

% \partial_t n = [2 . / ( (1 . +5 . *x . ~2) . * (1 + rho))-0.4]n 

% 

I 

clear all; 
close all; 

m=4000; "/onumber of points for space grid 

eps=.01; "/oepsilon 
rhoin=l.; "/oinitial rho 

dx=l/ni; X = (0:m)*dx; "/.space grid (0, l/m,2/m, . . . , 1) 
dt=25 . *dx~2/eps ; "/otime step 

dtsave = dt; Xinitialize array to save time steps 

tfinal = 15000 * dt; %final time 



weights = ones (numel (x) ) ; "/oweights for trapez rule integration 
weights (1) = 0.5; 
weights (end) = 0.5; 



27 



nd=exp(- (x- . 7) . ~2/eps) ; "/oinitial data 

nd=rhoin*nd/(dx*sum(nd) ) ; "/oinitial data normalized to have mass rhoin 

ndd=nd;R=x; "/oinitialize ndd and R 
count = 0; 

temps=0; 

savendd , : )=x(l : 10 :m) ; '/oinitialize arrary for saving nd at every 10th point in space 

while temps < tfinal 
count = count +1 
dt = min(dt, tfinal -temps); 
rho = dx * trapz(nd) 
savendC count , : )=nd(l : 10 :m) ; 
arrayt (count) = temps; 

arrayR(count) = ( (4-rho) / (5+5*rho) ) ' (1/2) ; 
R=2./((l.+5.*x."2) .*(1 + rho))-0.4; 

Rp=max(0,R)/eps; Rm=min(0,R)/eps; "/positive and negative part of R 

ndd = ((l+dt*Rp) .*nd ) . /(l-dt*Rm) ; '/.calculate n(t+dt) from n(t) 

nd=ndd ; 

temps = temps + dt; 

end 
f igure 

[X,Y] = meshgrid(arrayt,x(l:10:m)) ; 

[C,h]=contour(X,Y,savend' , [4 10 30] , 'LineColor' , 'black' ) ; 
text_handle = clabel(C,h); 
hold on; 

plot (arrayt ,arrayR, ' :k') ; 

xlabel('t') 

ylabel('x') 

figure ; 

plot (x,nd, 'k' ) ; 

xlabel('x') 

ylabel('n') 
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B Matlab code for Figure |5.2 



°/o solves equation 

% \partial_t n = [1 . / (1 • +x. "2)-dea- (al~2) . / (aa"2+x. "2)] n 

% 

% 

clear all; 
close all; 

m=4000; "/number of points for space grid 

eps=.01; 7oepsilon 
rhoin=l.; "/oinitial rho 

dx=l/m; X = (0:m)*dx; "/space grid (0, l/m,2/m, . . . , 1) 
dt=4500*dx~2/eps; "/time step 

dtsave = dt; "/initialize array to save time steps 

tfinal = 1000 * dt; "/final time 

weights = ones (numel (x) ) ; "/weights for trapez rule integration 
weights (1) = 0.5; 
weights (end) = 0.5; 

nd=exp(- (x- . 5) . "2/eps) ; "/initial data 

nd=rhoin*nd/(dx*sum(nd) ) ; "/initial data normalized to have mass rhoin 

ndd=nd;R=x; "/initialize ndd and R 

ndO=nd ; 
count = 0; 
temps=0; 

savendrenormd , : )=x(l : 10 :m) ; "/initialize arrary for saving nd at every 10th point in space 
"/ saveRd , : )=x(l : 10 :m) ; "/initialize arrary for saving R at every 10th point in space 

al=.55 
aa= . 5 
dea= . 245 

R=l./(l.+x.-2)-dea-(al-2) ./(aa~2+x. "2) ; 

while temps < tfinal 
count = count +1 

nd = exp(R*temps/eps) . *nd0; "/exponential formula to solve the ODE 
pp = nd/(dx * trapz(nd)); 
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Ik= dx * trapz (pp . *R) ; 

dt = minCdt, tfinal -temps); 
rho = dx * trapz (nd) 

savendrenormCcount , : )=iid(l : 10 :iii)/(dx * trapz(nd)); 
arr ay t (count) = temps; 

"/.calculating explicitly the position of the concentration point 

aq=(Ik+dea) 

bq=-l+al'2+(Ik+dea)*(l+aa'2) 
cq=-aa"2+al"2+ (Ik+dea) *aa"2 

xq2= (-bq- (bq"2-4*aq*cq) " (1/2) ) / (2*aq) 

xq = (xq2)-(l/2) 
arrayR( count) = xq; 

temps = temps + dt; 



end 

figure 

[X,Y] = meshgrid(arrayt ,x(l : 10 :m) ) ; 

[C,h] =contour (X, Y, savendrenorm' , [2 8 16] , 'LineColor' , 'black') ; 
text_handle = clabel(C,h); 
hold on; 

plot (arrayt , arrayR, ' : k ' ) ; 

xlabel('t') 
ylabel('x') 

f igure ; 

plot (x,nd, 'k' ) ; 

xlabel('x') 

ylabel('n') 
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